#include "Astronomy/Nutation.hpp"
#include "Astronomy/Transformation.hpp"



const std::array<HighNutationCoefficient, 77> CNutation::HighNutationTable{
{
        //l  l' F  D  Om         A        A'      A"        B      B'     B"
       {  0, 0, 0, 0, 1, -172064161, -174666,  33386, 92052331,  9086, 15377},
       {  0, 0, 2,-2, 2,  -13170906,   -1675, -13696,  5730336, -3015, -4587},
       {  0, 0, 2, 0, 2,   -2276413,    -234,   2796,   978459,  -485,  1374},
       {  0, 0, 0, 0, 2,    2074554,     207,   -698,  -897492,   470,  -291},
       {  0, 1, 0, 0, 0,    1475877,   -3633,  11817,    73871,  -184, -1924},
       {  0, 1, 2,-2, 2,    -516821,    1226,   -524,   224386,  -677,  -174},
       {  1, 0, 0, 0, 0,     711159,      73,   -872,    -6750,     0,   358},
       {  0, 0, 2, 0, 1,    -387298,    -367,    380,   200728,    18,   318},
       {  1, 0, 2, 0, 2,    -301461,     -36,    816,   129025,   -63,   367},
       {  0,-1, 2,-2, 2,     215829,    -494,    111,   -95929,   299,   132},

       {  0, 0, 2,-2, 1,     128227,     137,    181,   -68982,    -9,    39},
       { -1, 0, 2, 0, 2,     123457,      11,     19,   -53311,    32,    -4},
       { -1, 0, 0, 2, 0,     156994,      10,   -168,    -1235,     0,    82},
       {  1, 0, 0, 0, 1,      63110,      63,     27,   -33228,     0,    -9},
       { -1, 0, 0, 0, 1,     -57976,     -63,   -189,    31429,     0,   -75},
       { -1, 0, 2, 2, 2,     -59641,     -11,    149,    25543,   -11,    66},
       {  1, 0, 2, 0, 1,     -51613,     -42,    129,    26366,     0,    78},
       { -2, 0, 2, 0, 1,      45893,      50,     31,   -24236,   -10,    20},
       {  0, 0, 0, 2, 0,      63384,      11,   -150,    -1220,     0,    29},
       {  0, 0, 2, 2, 2,     -38571,      -1,    158,    16452,   -11,    68},

       {  0,-2, 2,-2, 2,      32481,       0,      0,   -13870,     0,     0},
       { -2, 0, 0, 2, 0,     -47722,       0,    -18,      477,     0,   -25},
       {  2, 0, 2, 0, 2,     -31046,      -1,    131,    13238,   -11,    59},
       {  1, 0, 2,-2, 2,      28593,       0,     -1,   -12338,    10,    -3},
       { -1, 0, 2, 0, 1,      20441,      21,     10,   -10758,     0,    -3},
       {  2, 0, 0, 0, 0,      29243,       0,    -74,     -609,     0,    13},
       {  0, 0, 2, 0, 0,      25887,       0,    -66,     -550,     0,    11},
       {  0, 1, 0, 0, 1,     -14053,     -25,     79,     8551,    -2,   -45},
       { -1, 0, 0, 2, 1,      15164,      10,     11,    -8001,     0,    -1},
       {  0, 2, 2,-2, 2,     -15794,      72,    -16,     6850,   -42,    -5},

       {  0, 0,-2, 2, 0,      21783,       0,     13,     -167,     0,    13},
       {  1, 0, 0,-2, 1,     -12873,     -10,    -37,     6953,     0,   -14},
       {  0,-1, 0, 0, 1,     -12654,      11,     63,     6415,     0,    26},
       { -1, 0, 2, 2, 1,     -10204,       0,     25,     5222,     0,    15},
       {  0, 2, 0, 0, 0,      16707,     -85,    -10,      168,    -1,    10},
       {  1, 0, 2, 2, 2,      -7691,       0,     44,     3268,     0,    19},
       { -2, 0, 2, 0, 0,     -11024,       0,    -14,      104,     0,     2},
       {  0, 1, 2, 0, 2,       7566,     -21,    -11,    -3250,     0,    -5},
       {  0, 0, 2, 2, 1,      -6637,     -11,     25,     3353,     0,    14},
       {  0,-1, 2, 0, 2,      -7141,      21,      8,     3070,     0,     4},

       {  0, 0, 0, 2, 1,      -6302,     -11,      2,     3272,     0,     4},
       {  1, 0, 2,-2, 1,       5800,      10,      2,    -3045,     0,    -1},
       {  2, 0, 2,-2, 2,       6443,       0,     -7,    -2768,     0,    -4},
       { -2, 0, 0, 2, 1,      -5774,     -11,    -15,     3041,     0,    -5},
       {  2, 0, 2, 0, 1,      -5350,       0,     21,     2695,     0,    12},
       {  0,-1, 2,-2, 1,      -4752,     -11,     -3,     2719,     0,    -3},
       {  0, 0, 0,-2, 1,      -4940,     -11,    -21,     2720,     0,    -9},
       { -1,-1, 0, 2, 0,       7350,       0,     -8,      -51,     0,     4},
       {  2, 0, 0,-2, 1,       4065,       0,      6,    -2206,     0,     1},
       {  1, 0, 0, 2, 0,       6579,       0,    -24,     -199,     0,     2},

       {  0, 1, 2,-2, 1,       3579,       0,      5,    -1900,     0,     1},
       {  1,-1, 0, 0, 0,       4725,       0,     -6,      -41,     0,     3},
       { -2, 0, 2, 0, 2,      -3075,       0,     -2,     1313,     0,    -1},
       {  3, 0, 2, 0, 2,      -2904,       0,     15,     1233,     0,     7},
       {  0,-1, 0, 2, 0,       4348,       0,    -10,      -81,     0,     2},
       {  1,-1, 2, 0, 2,      -2878,       0,      8,     1232,     0,     4},
       {  0, 0, 0, 1, 0,      -4230,       0,      5,      -20,     0,    -2},
       { -1,-1, 2, 2, 2,      -2819,       0,      7,     1207,     0,     3},
       { -1, 0, 2, 0, 0,      -4056,       0,      5,       40,     0,    -2},
       {  0,-1, 2, 2, 2,      -2647,       0,     11,     1129,     0,     5},

       { -2, 0, 0, 0, 1,      -2294,       0,    -10,     1266,     0,    -4},
       {  1, 1, 2, 0, 2,       2481,       0,     -7,    -1062,     0,    -3},
       {  2, 0, 0, 0, 1,       2179,       0,     -2,    -1129,     0,    -2},
       { -1, 1, 0, 1, 0,       3276,       0,      1,       -9,     0,     0},
       {  1, 1, 0, 0, 0,      -3389,       0,      5,       35,     0,    -2},
       {  1, 0, 2, 0, 0,       3339,       0,    -13,     -107,     0,     1},
       { -1, 0, 2,-2, 1,      -1987,       0,     -6,     1073,     0,    -2},
       {  1, 0, 0, 0, 2,      -1981,       0,      0,      854,     0,     0},
       { -1, 0, 0, 1, 0,       4026,       0,   -353,     -553,     0,  -139},
       {  0, 0, 2, 1, 2,       1660,       0,     -5,     -710,     0,    -2},

       { -1, 0, 2, 4, 2,      -1521,       0,      9,      647,     0,     4},
       { -1, 1, 0, 1, 1,       1314,       0,      0,     -700,     0,     0},
       {  0,-2, 2,-2, 1,      -1283,       0,      0,      672,     0,     0},
       {  1, 0, 2, 2, 1,      -1331,       0,      8,      663,     0,     4},
       { -2, 0, 2, 2, 2,       1383,       0,     -2,     -594,     0,    -2},
       { -1, 0, 0, 0, 2,       1405,       0,      4,     -610,     0,     2},
       {  1, 1, 2,-2, 2,       1290,       0,      0,     -556,     0,     0}
}};
//中精度章动计算表;
const std::array<MediumNutationCoefficient, 10> CNutation::MediumNutationTable{
{
        {  2.1824,  -33.75705,  36e-6, -1720, 920},
        {  3.5069, 1256.66393,  11e-6,  -132,  57},
        {  1.3375,16799.41820, -51e-6,   -23,  10},
        {  4.3649,  -67.51410,  72e-6,    21,  -9},
        {    0.04, -628.30200,      0,   -14,   0},
        {    2.36, 8328.69100,      0,     7,   0},
        {    3.46, 1884.96600,      0,    -5,   2},
        {    5.44,16833.17500,      0,    -4,   2},
        {    3.69,25128.11000,      0,    -3,   0},
        {    3.55,  628.36200,      0,     2,   0}

 }};




NutationAngles CNutation::HighPrecision(double T, const int NumberOfCycleDays) noexcept
{
    const double Tsquared{ T * T };
    const double Tcubed{ Tsquared * T };
    const double T4{ Tcubed * T };

    const double M = 485868.249036 + 1717915923.2178 * T + 31.8792 * Tsquared + 0.051635 * Tcubed - 0.00024470 * T4;
    const double Mprime = 1287104.79305 + 129596581.04810 * T - 0.55320 * Tsquared - 0.000136 * Tcubed - 0.00001149 * T4;
    const double F = 335779.526232 + 1739527262.8478 * T - 12.7512 * Tsquared - 0.001037 * Tcubed + 0.00000417 * T4;
    const double D = 1072260.70369 + 1602961601.2090 * T - 6.37060 * Tsquared + 0.006593 * Tcubed - 0.00003169 * T4;
    const double Omega = 450160.398036 - 6962890.5431000 * T + 7.47220 * Tsquared + 0.007702 * Tcubed - 0.00005939 * T4;

    double dL = 0, dE = 0;

    //周期项取和计算;
    for (const auto& coeff : CNutation::HighNutationTable)
    {
        double sita = coeff.M * M + coeff.Mprime * Mprime + coeff.F * F + coeff.D * D + coeff.omega * Omega;
        sita = Transformation::ArcsecToRadians(sita);

        
        if (NumberOfCycleDays != 0)
        {
            const double CurrentSita = (1717915923.2178 * coeff.M + 129596581.0481 * coeff.Mprime + 1739527262.8478 * coeff.F + 1602961601.2090 * coeff.D + 6962890.5431 * coeff.omega);
            const double CurrentArcsec = Transformation::RadiansToArcsec(36526.0 * 2.0 * Constants::PI);
            const double Current = CurrentArcsec / CurrentSita;

            //只算周期大于NumberOfCycleDays天的项;
            if (Current < NumberOfCycleDays)
                continue;
        }

        const double sinArg = sin(sita);
        const double cosArg = cos(sita);

        dL += (coeff.sincoeff1 + coeff.sincoeff2 * T) * sinArg + coeff.sincoeff3 * cosArg;
        dE += (coeff.coscoeff1 + coeff.coscoeff2 * T) * cosArg + coeff.coscoeff3 * sinArg;
    }

    // return{ dL /= 10000000 * rad, dE /= 10000000 * rad };//返回IAU2000B章动值, dL是黄经章动,dE是交角章动
    // 
     //章动表的系数10000000.0;
    return { Transformation::ArcsecToRadians(dL / 10000000.0), Transformation::ArcsecToRadians(dE / 10000000.0) };
}

NutationAngles CNutation::MediumPrecision(double T) noexcept
{
     const double Tsquared{ T * T };
 
     double value = -1.742 * T;

     double dL = 0, dE = 0;
     for (auto& coffe : CNutation::MediumNutationTable)
     {
         double sita = coffe.T + coffe.T1 * T + coffe.T2 * Tsquared;
         dL += (coffe.sincoeff1 + value) * sin(sita);
         dE += coffe.coscoeff1 * cos(sita);

         value = 0;
     }

     // 章动表的系数100.0;
     return { Transformation::ArcsecToRadians(dL / 100.0), Transformation::ArcsecToRadians(dE / 100.0) };
}
double CNutation::MediumPrecisionLongitude(double T) noexcept
{
    const double Tsquared{ T * T };

    double value = -1.742 * T;

    double dL = 0;
    for (auto& coffe : CNutation::MediumNutationTable)
    {
        double sita = coffe.T + coffe.T1 * T + coffe.T2 * Tsquared;
        dL += (coffe.sincoeff1 + value) * sin(sita);
        value = 0;
    }

    // 章动表的系数100.0;
    return  Transformation::ArcsecToRadians(dL / 100.0);
}

Coordinate3D CNutation::EquatorNutation(const Coordinate3D& equator, const double epsilon, const NutationAngles& nutation) noexcept
{
    Coordinate3D result{ equator };

    const double sinE = sin(epsilon);
    const double cosX = cos(equator.X);
    const double sinX = sin(equator.X);

    result.X += (cos(epsilon) + sinE * sinX * tan(equator.Y)) * nutation.longitude - cosX * tan(equator.Y) * nutation.obliquity;
    result.Y += sinE * cosX * nutation.longitude + sinX * nutation.obliquity;

    result.X = Transformation::MapTo0To2PIRange(result.X);
    return result;
}